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ABSTRACT The glycocalyx is a sugar-rich layer located at the luminal part of the endothelial cells. It is involved in key meta¬ 
bolic processes and its malfunction is related to several diseases. To understand the function of the glycocalyx, a molecular level 
characterization is necessary. In this article, we present large-scale molecular-dynamics simulations that provide a comprehen¬ 
sive description of the structure and dynamics of the glycocalyx. We introduce the most detailed, to-date, all-atom glycocalyx 
model, composed of lipid bilayer, proteoglycan dimers, and heparan sulfate chains with realistic sequences. Our results reveal 
the folding of proteoglycan ectodomain and the extended conformation of heparan sulfate chains. Furthermore, we study the 
glycocalyx response under shear flow and its role as a flypaper for binding fibroblast growth factors (FGFs), which are involved 
in diverse functions related to cellular differentiation, including angiogenesis, morphogenesis, and wound healing. The simula¬ 
tions show that the glycocalyx increases the effective concentration of FGFs, leading to FGF oligomerization, and acts as a lever 
to transfer mechanical stimulus into the cytoplasmic side of endothelial cells. 


INTRODUCTION 

The functioning of the human metabolism requires transfer 
of nutrients, gas molecules, and hormones from the blood¬ 
stream into the target organs. Such exchange between cir¬ 
culating blood and the tissues occurs across a thin layer of 
endothelial cells located in the luminal surface of the blood 
vessels, known as the endothelium. The endothelium forms 
the inner cellular lining covering the entire vascular 
network, from the wide heart arteries to the smallest capil¬ 
laries (1). Far from being an inert tissue, the endothelium 
is very versatile, and it is conveniently heterogeneous (2), 
adapting itself to the local environment. Moreover, the 
endothelium participates actively in critical biological 
processes such as blood cell trafficking, angiogenesis, and 
immune response ( 3 ), and its malfunction is related to 
several diseases (4,5). 

The glycocalyx is the outermost part of the endothelium. 
Located at the surface of endothelial cells, the glycocalyx is 
the first vascular barrier in direct contact with blood. Fig. 1 
shows an atomic model of the glycocalyx in the context of 
the circulatory system. As its name states, the glycocalyx 
is a sugar-rich layer, with carbohydrates associated by cova¬ 
lent links forming oligosaccharide chains that can be either 
long and linear or short and branched. The presence of 
anionic oligosaccharides, such as heparan sulfate, provides 
the glycocalyx layer with a characteristic negative charge. 
Most oligosaccharides are covalently attached to membrane 
proteins, but it is also possible to find that some are free and 
remain adsorbed in the glycocalyx region by noncovalent 
interactions. Furthermore, the glycocalyx can retain plasma 
proteins recruited from extracellular space and its thickness 
can be lessened by enzymatic and shear-induced shedding, 
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producing a dynamic equilibrium between its composition 
and the flowing blood (6). Due to the variety of its com¬ 
position and its exposed location, the glycocalyx plays 
two critical functions for the endothelium: as an interface 
between chemical signals in the blood and biochemical 
pathways inside cells; and as a barrier against chemical or 
mechanical damage. 

In early studies of the endothelium, the glycocalyx had 
remained unnoticed because it was difficult to visualize: 
the tiny sugar coat was washed out by the preparation 
procedures involved in conventional light and electron mi¬ 
croscopy. The first visual evidences of the glycocalyx 
came through the employment of cationic dyes, such as 
ruthenium red, which were strongly adsorbed by the nega¬ 
tively charged glycocalyx (7). Further improvements in 
staining protocols as well as in spectroscopic methods 
have allowed obtaining static images of the glycocalyx 
(1). Since then, research in glycocalyx has flourished and 
now it is investigated as a therapeutic target against several 
pathologies, such as atherosclerosis (8) and diabetes (9). 
Furthermore, recent experimental studies hypothesize a 
long-range effect of glycocalyx on water ordering (10), 
which may mediate cell surface interactions. 

Due to the experimental difficulty of tracing the elusive 
atomic events inside the glycocalyx, researchers have turned 
to computer simulations (11-16). Those models were 
mostly based on continuum or coarse-grained approaches 
but refined to include the peculiarities of the glycocalyx, 
such as the negatively charged surface or the steric resis¬ 
tance of oligosaccharide chains. Although computer simula¬ 
tions have contributed to contemporary understanding of the 
glycocalyx, structural information is still much needed to 
advance the field. We turn to all-atom molecular dynamics 
(MD) simulations, which provide the most detailed and dy¬ 
namic representation of biomolecular systems (17,18). 
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FIGURE 1 Glycocalyx in the circulatory system, (a) A cartoon represen¬ 
tation of a blood vessel: ( orange circles ) endothelial cells; ( green brush) 
glycocalyx. ( b ) All-atom model of the glycocalyx. {Green licorice repre¬ 
sentation) Sugar chains. {Red surfaces) Transmembrane protein; {van der 
Waals beads) lipid bilayer. Water and ions are not shown. 


The second group relates to glycocalyx simulations under shear flow 
(Sim33-35). 

The third group lists simulations used to study binding of proteins to the 
glycocalyx (Sim36^f0). 


Building the glycocalyx 

We chose Syndecan-4 (Syn-4) proteoglycan and heparan sulfate (HS) sugar 
chains for modeling the glycocalyx as described in the Results and Discus¬ 
sion. Below we provide information about assembly protocol, simulation 
time, and system size. 

We initially divided the glycocalyx structure into three parts: Syndecan-4 
ectodomain (Syn-4 Ect) linked with sugar chains; Syndecan-4 transmem¬ 
brane (Syn-4 Tmb) dimer embedded into a lipid bilayer; and Syndecan-4 
cytoplasmic (Syn-4 Cyt) dimer. Fig. SI in the Supporting Material shows 
the three separate systems. After reaching equilibration for each part, the 
sections were glued together into a single glycocalyx system and equili¬ 
brated further. 

For Syn-4 Ect simulations (Sim 1-24), our starting structure was a fully 
stretched ectodomain joined with three extended 100-residue sugar chains. 
Subsequently, the structure was annealed as follows: the system was heated 
at 1000 K for 10 ns, then cooled down to 310 K at a rate of —100 K/l ns. 
After that, the structure was checked for stereochemical correctness (19). 
The final annealing stage was an equilibration of 100 ns, summing a total 
of 116 ns per annealing cycle. The tested range of temperatures does not 
allow employing all-atom water models; therefore, an implicit solvent 
model was used. We simulated four different HS sequences (see Fig. 2 
and the Results and Discussion). Six different conformations were tested 
for each HS sequence. The chains orientations were random, but avoided 
overlapping contacts between residues. All together, the simulation time 
spent in annealing cycles was 

(4 HS options) x (6 conformations) 
x (116 ns per cycle) = 2.784^8. 


Here, we report large-scale, all-atom MD simulations of 
the glycocalyx layer. To begin, we describe a procedure to 
build all-atom models of glycocalyx. We present structural 
insight obtained while building the model: the folding of 
the ectodomain of proteoglycans, the extended conforma¬ 
tion of sugar chains, and the negligible effect of sulfonation 
charge on the chain length. Finally, our model was simu¬ 
lated under nonequilibrium conditions to address two 
fundamental questions about the glycocalyx: its response 
to mechanical stimuli under shear flow, and its ability to 
retain fibroblast growth factor (FGF) proteins. 

METHODS 

MD simulations compute the time-evolution of a many-body system by solv¬ 
ing Newton’s equations of motion. Two compulsory requirements for MD 
simulations are the initial atomic positions of particles, and the force field 
describing the interactions between particles. The glycocalyx is intrinsically 
disordered and it is not possible to obtain a minute atomic structure by NMR 
or crystallography; thus, the atomic models relevant to this discussion were 
built based on the available structural information as detailed below. 

The MD simulations are organized into three groups: 

The first group involves simulations (Sim) needed to build the glycocalyx 
structure (Siml-32 in Table 1): i.e., annealing of ectodomain (Siml-24), 
equilibration of transmembrane (Sim25-26) and cytoplasmic (Sim27-28) 
domains, and the assembly of 5,700,000-atom glycocalyx systems 
(Sim29-32). 


For Syn-4 Tmb simulations, we mapped two Syn-4 backbones into a re¬ 
ported Glycophorin-A dimer structure (20) (PDB:lAFO; see Results and 
Discussion). Weinbaum et al. (15) have proposed a hexagonal symmetry 
for glycocalyx. In accordance with such model, Syn-4 transmembrane di¬ 
mers were embedded into a POPC lipid bilayer of hexagonal periodicity 
and 820 nm 2 area. The lipid membranes were generated from a preequili¬ 
brated membrane patch provided with the MEMBRANE BUILDER plugin 
available in the software package VMD (21). Two systems were built, the 
first containing one Syn-Tmb dimer, the second with three Syn-Tmb dimers 
equally spaced among themselves and their periodic images. Then, the sys¬ 
tems were covered with water boxes resulting in 11-nm height each. Both 
systems were equilibrated in NPT ensemble for 5 ns (Sim25-26). 

For Syn-4 Cyt simulations, we solvated and ionized a Syn-4 cytoplasmic 
dimer (22) (PDB:1EJP) using a water box of 8 x 8 x 9 nm 3 . Then, the 
system was equilibrated for 5 ns in NPT ensemble (Sim27), followed by 
5 ns in NVT (Sim28). 

The final stage is to join all three parts: the transmembrane domain con¬ 
tains residues 141-171; the ecto domain, 1-141; and the cytoplasmic domain, 
171-198 (see Fig. 2). The repeating residues in the transmembrane part were 
used as guiding positions to attach the other domains. Special care was taken 
to avoid sugar rings either intersecting each other across the periodic images 
or overlapping the membrane region. Once the three parts were joined, topol¬ 
ogy was reconstructed around the new bonds and the entire structure was 
checked for chirality correctness. Two systems were built, with one and three 
Syn-4 dimers, then solvated and ionized to 100 mM NaCl concentration. The 
final systems are composed of 5,700,000 atoms, with 820 nm 2 area and 69 nm 
height. MD equilibrations proceeded as follows: 3.5 ns in NPT ensemble 
(Sim29-30) and 5 ns in NVT (Sim31-32). The final equilibrated structures 
for both glycocalyx systems are provided in the Supporting Material. 
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TABLE 1 Summary of MD simulations (Sim) 


Label 

Type' 1 

Atoms (x 1000) 

Time (ns) 

Solvent model 

Starting structure 1 

Siml-6 

Annealing 

9 

116 

I 

Syn-4 Ect and [NA-NA/NS-NS] sequence 

Sim7—12 

Annealing 

9 

116 

I 

Syn-4 Ect and NA sequence 

Siml3-18 

Annealing 

9 

116 

I 

Syn-4 Ect and NS sequence 

Siml9-24 

Annealing 

9 

116 

I 

Syn-4 Ect and NAH sequence 

Sim25 

Eq-NPT 

955 

5 

E 

Syn-4 Tmb dimer 

Sim26 

Eq-NPT 

955 

5 

E 

3 Syn-4 Tmb dimers 

Sim27 

Eq-NPT 

58 

5 

E 

Syn-4 Cyt dimer 

Sim28 

Eq-NVT 

58 

5 

E 

Sim27 

Sim29 

Eq-NPT 

5700 

3.5 

E 

Gccx (Glycocalyx) (1 Syn-4 dimer) 

Sim30 

Eq-NPT 

5700 

3.5 

E 

Gccx (Glycocalyx) (3 Syn-4 dimers) 

Sim31 

Eq-NVT 

5700 

5 

E 

Sim29 

Sim32 

Eq-NVT 

5700 

5 

E 

Sim30 

Sim33 

Flow-NVT 

5950 

0.6 

E 

Sim31 

Sim34 

Flow-NVT 

5950 

1.5 

E 

Sim31 

Sim35 

Flow-NVT 

5950 

10 

E 

Sim31 

Sim36 

NVT 

20 

10 

I 

10 FGF-2 

Sim37 

GSMD-NVT 

20 

10 

I 

Sim36 and 1 Syn-4 dimer grid 

Sim38 

GSMD-NVT 

20 

10 

I 

Sim36 and 3 Syn-4 dimer grid 

Sim39 

NVT 

20 

10 

I 

Sim37 

Sim40 

NVT 

20 

10 

I 

Sim38 


Table presents the simulation effort for building (Sim 1-32), flow (Sim33-35), and binding simulations (Sim36-40). 
a Eq, equilibration; GSMD, Grid-SMD; NVT, constant volume; NPT, constant pressure. 
b 7, generalized Bom implicit solvent; E, TIP3P explicit solvent. 
c Syn, Syndecan; Ect, ectodomain; Gccx, glycocalyx. 


Shear-stress effect on glycocalyx conformation 

We study the effect of flow on the glycocalyx conformation by applying 
forces on the water molecules to mimic shear-stress. Our starting structure 
was the last frame of simulation Sim32. It is worth noting that the system is 
periodic in the three axes. Therefore, shear-stress applied on the ecto¬ 
plasmic side would propagate to the cytoplasmic side through the z-periodic 
boundary and not across the Syn-4 dimers embedded in the lipid membrane. 
To avoid this caveat, the system was modified along the z axis by adding a 
layer of eight graphene sheets, located 51 nm above the lipid membrane. 
The graphene layer has a thickness of 2 nm, the same area as the lipid mem¬ 
brane, and it is periodic in x- and y axes, creating an effective division 
between two compartments: the ectoplasmic compartment from the top 
of lipid membrane to the bottom of the graphene layer; and the cytoplasmic 
compartment from the bottom of lipid membrane to the top of the graphene 
layer. All HS sugar chains are contained in the ectoplasmic compartment. 


The graphene layers were created using the INORGANIC BUILDER plu¬ 
gin (23) in the software VMD. The system containing the glycocalyx and 
graphene barrier was minimized for 1000 steps and equilibrated for 
0.2 ns in NPT ensemble. Graphene atoms were restrained in the x- and y 
axes to keep the area constant. For production runs (Sim33-35), the gra¬ 
phene layer remains fixed; therefore, shear-stress cannot cross the z-peri- 
odic boundary. The forces are applied on the water oxygens located in 
the upper half of the ectoplasmic compartment. The direction of the force 
is along the x axis, parallel to the lipid membrane. 

Glycocalyx binding 

We study the effects of the glycocalyx on flowing proteins using grid- 
steered molecular dynamics (Grid-SMD) simulations. In Grid-SMD, the 
key feature of the glycocalyx, namely its electrostatic potential, is mapped 
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FIGURE 2 Sequences for transmembrane pro¬ 
tein and heparan sulfate chains, (a) Amino-acid 
sequence of Syndecan-4. Residues located in the 
transmembrane region are underlined. Residues 
from 1 to 141 are located in the extracellular re¬ 
gion. Proline residues ( highlighted) contribute to 
disrupt secondary structure formation in the ecto¬ 
domain. Three serine residues ( italicized , high¬ 
lighted with arrows ), are attaching points for 
heparan sulfate (HS) chains. ( b-e ) Four sugar 
sequences of HS chains (see text for details). ( b ) 
Intercalation of sulfonated and nonsulfonated HS. 
(c) Nonsulfonated HS. ( d) Fully sulfonated HS. 
( e ) Protonated HS. ( Symbols ) Different monosac¬ 
charides. 
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into a three-dimensional grid. Overall, the protocol for binding simulations 
is similar to the molecular-dynamics flexible fitting technique (24), where 
biomolecules are force-fitted into cryo-electron microscopy grids. In our 
case, the flowing proteins are force-driven into glycocalyx grids. The grids 
were calculated using the PMEPOT (25) plugin available with the software 
VMD. PMEPOT calculates the electrostatic potential from a MD trajectory 
by recomputing the long-range electrostatic forces via particle-mesh Ewald 
(PME) (26). Two grids were computed, one for the system with one Syn-4 
dimer, the other for the system with three Syn-4 dimers. We characterized 
the glycocalyx grids by the average electrostatic potential of 5-ns MD 
equilibration trajectories (Sim31 and 32), using only the atomic charges 
of HS chains, whereas other components, such as Syn-4 dimers, lipid mem¬ 
brane, and electrolytic solution were not included. This simplification was 
necessary to provide a smooth potential that drives flowing proteins near 
highly charged sugar locations within dozens of nanoseconds. In the Sup¬ 
porting Material, we include references to a Grid-SMD tutorial and the gly¬ 
cocalyx grids. 

The flowing protein is fibroblast growth factor-2 (FGF-2) (PDB:1BFB). 
We refer to the Results and Discussion for more details about FGF-2. To 
induce flow, a force of 1 pN was applied to the a-carbons of FGF-2, with 
the direction of the force being along the x axis. For binding simulations, 
we performed a three-step MD cycle. 

For the first part (Sim36), 10 all-atom FGF-2 proteins were randomly 
distributed on a region from 14 to 39 nm above the membrane location. 
The minimum distance between FGF-2’s centers of mass was 5 nm. Peri¬ 
odic boundary conditions were assumed, and the periodic lengths matched 
the glycocalyx grid dimensions. Then, FGF-2 molecules were allowed to 
flow for 10 ns. 

For the second part, two separate 10-ns MD simulations were performed, 
one with the single Syn-4 dimer glycocalyx grid (Sim37) and the other with 
the three Syn-4 dimer glycocalyx grid (Sim38). The starting point was the 
last frame of Sim36, but now the glycocalyx grids were active, affecting 
FGF-2 flow dynamics. 

For the third part, both Sim37 and 38 simulations were continued for 
10 ns, but the glycocalyx grids were turned off (Sim39 and 40, respec¬ 
tively). For all three parts, implicit solvent was used. To maintain FGF-2 
structure, secondary structure restraints were used and the steering forces 
derived from the grid were scaled up by a factor of 0.05. The system size 
was 20,180 atoms, consisting of only 10 FGF-2 molecules. 


6,000,000 atom systems, we used the memory-optimized version of 
NAMD (35). 

The analysis was performed using VMD (Visual Molecular Dynamics, 
University of Illinois at Urbana-Champaign, Champaign, IL), MATLAB 
(The Math Works, Natick, MA), and PYTHON (Python Software Founda¬ 
tion, Wilmington, DE) scripts. All MD snapshots were made with VMD. 
The secondary structure plots (Fig. 3 a), were computed using the TimeLine 
plugin in VMD. To calculate the velocity profiles, the system was divided in 
layers of 1 nm height along the z axis; for each layer, the velocities in the x 
axis of water oxygens were averaged. The shear stresses (r) were computed 
from the resulting velocity profiles as 


<5v 

T = M &’ 


where p is the viscosity of the TIP3P water model (36) (0.321 mPa« s); and 
5v/5z is the gradient of the x velocity along the z axis. The persistence length 
(Pc) was calculated using the relationship 

( e 2e 2 ) = 2P C L C - 2P 2 c [1 - exp (~L C /P C )\, 


where ( ele 2 ) is the averaged squared end-to-end distance, L c is the contour 
length, and Pc is the persistence length (37). 
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Protocol details 

All MD simulations were performed using the software program NAMD 
2.7 (27). We used the CHARMM biomolecular force field (28) and the 
TIP3P water model (29). Structures were generated using the PSFgen utility 
provided with the VMD program. The latest additions to CHARMM 
include parameters for sugars (30-32).To assure compatibility with 
PSFgen, atom types were renamed from the six-letter code given in the 
CHARMM force field into a four-letter code. Sulfamate group parameters 
(R-NHSO _1 3 ) needed for sulfonated GlcNAc were taken from Huige and 
Altona (33). 

For MD simulations with implicit solvent, namely annealing (Sim 1-24) 
and binding (Sim36-37) simulations, we used 1-fs timestep and flexible 
bonds. We employed the generalized Bom implicit solvent (34) that accu¬ 
rately represents the dielectric screening between solute molecules in elec¬ 
trolytic solution. The van der Waals and electrostatic interactions were 
calculated using a cutoff of 16 A with a switching function starting at 
15 A and 0.1 M electrostatic screening. Temperature was controlled using 
a Langevin thermostat. For MD simulations with explicit water (Sim25- 
35), we used a 2-fs timestep, rigid bonds, and PME electrostatics with a 
grid density of 1/A 3 . The van der Waals interactions were calculated using 
a cutoff of 12 A with a switching function starting at 10 A. Simulations in 
NPT ensemble were kept at 1 atm pressure and 310 K using a Langevin ther¬ 
mostat and a Nose-Hoover Langevin piston. Simulations in the NVT 
ensemble were kept at 310 K using the Andersen thermostat. To simulate 


C 



FIGURE 3 Ectodomain folding, (a) Plot shows the time evolution of sec¬ 
ondary stmcture for Syn-4 ectodomain during the annealing cycle. Abscissa 
and ordinate axes show time and residue number, respectively. (Vertical 
lines) The three annealing stages: 10 ns heating at 1000 K, 6 ns cooling 
from 1000 to 310 K, and final 100-ns equilibration at 310 K. (Shaded 
bands) Appearance of a-helical stmcture. (b) RMSD of a:-carbons of 
Syn-4 during equilibration at 310 K. (c) (Snapshot) Last frame of simula¬ 
tion, with (shaded ribbons) a-helices, (solid spheres) proline residues, 
and (lines) HS chains. 
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RESULTS AND DISCUSSION 
Building glycocalyx system 

The glycocalyx layer is composed of sugar chains cova¬ 
lently linked to transmembrane proteins located at the sur¬ 
face of endothelial cells. A complete all-atom model of 
the glycocalyx is notoriously difficult to build: the oligosac¬ 
charide chains are very heterogeneous in structure and 
sequence to be resolved by crystallization or NMR and 
only low-resolution images are available (38-40). To over¬ 
come this problem, most glycocalyx models have only 
considered some key features such as negative charge or 
polymer hindrance (11-16), while disregarding other 
equally important features such as nonuniform charge distri¬ 
bution, realistic sugar chain length, or interactions among 
the different biomolecular components. 

We aimed to build the most detailed all-atom glycocalyx 
model using structural information. Some degree of simpli¬ 
fication is required for a system with such high structural 
heterogeneity, but we did preserve the molecular compo¬ 
nents that are most prevalent in glycocalyx. The most abun¬ 
dant oligosaccharides in the glycocalyx are heparan sulfate 
(50-90%) and chondroitin sulfate (10-20%) (6). Accord¬ 
ingly, our glycocalyx considers heparan sulfate (HS) chains. 
The anchoring transmembrane proteins for HS chains are 
syndecans. We chose Syndecan-4 (Syn-4), which is widely 
distributed among all endothelial cells, in contrast to other 
syndecans that are localized in specific tissues (41). Syn-4 
is a proteoglycan of 198-amino-acid length, containing 
extracellular, transmembrane, and cytoplasmic domains. 
The extracellular domain (i.e., the ectodomain) contains 
three points for attaching HS molecules (42). Fig. 2 presents 
the employed Syn-4 and HS sequences. 

Based on data from electron microscopy images (38), 
Weinbaum et al. (15) have presented a quasiperiodic hexag¬ 
onal arrangement for a glycocalyx with a spacing of 100 nm 
between proteoglycan clusters. Also, HS chains are long 
linear polymers of -100-200 sugar residues, corresponding 
to 50-100-nm length (1,43). An all-atom model with 
such spacing and dimensions would demand -100,000,000 
atoms, which is at the limit of contemporary MD capabil¬ 
ities. To reduce the computational cost of building the struc¬ 
ture, the system was initially split into three smaller pieces: 
Syn-4 ectodomain linked with HS chains; Syn-4 transmem¬ 
brane dimer embedded into a lipid bilayer; and Syn-4 cyto¬ 
plasmic dimer (see Methods). 

The first topology block is Syn-4 ectodomain joined with 
three HS chains. The HS length is assumed to be 100 sugar 
residues (44). The HS chains are covalently attached to three 
serine residues at positions 39, 61, and 63 in the ectodomain. 
The linker between HS and ectodomain has been described 
in other publications (6,45). In short, HS and Syn-4 are 
joined through a tetrasaccharide GlcA-[/j I -3]-Gal-[/31 -3]- 
Gal-[/?l^l]-Xyl (see Fig. 2). At the reducing end ( right- 
hand side of Fig. 2, b-e), a xylose sugar attaches to a serine 


amino acid through (3 1 glycosidic linkage. The HS chain 
grows at the nonreducing end ( left-hand side of Fig. 2, 
b-e) by adding either GlcA-[/31-4]-Glc-[al-4] or IdoA- 
[a I -4]-GIcNAc-[ a I -4 ] disaccharide units. The former 
disaccharide appears in nonsulfonated form, also known 
as the NA region. The latter disaccharide can be either 
partially sulfonated (NA/NS region) or highly sulfonated 
(NS region). A HS chain usually contains alternating NA, 
NA/NS, and NS regions of variable length, and the specific 
sequence seems to correlate with the cell type origin (45). 
Although increasing evidence suggests precise sulfonation 
patterns with strong specificity for particular proteins (46), 
current experimental techniques can only partially resolve 
HS sequences (47). To provide the closest model to a real 
glycocalyx, the sulfonation patterns for NS and NA/NS re¬ 
gions were chosen to bind fibroblast growth factor and 
antithrombin, respectively (45) (see Fig. 2 b). At the non¬ 
reducing end, we adopted the sulfonation pattern proposed 
by Wu and Lech (48). 

Due to the presence of 12 proline residues and a lack of 
homology structures, the ectodomain is usually depicted 
in a centipedelike representation; that is, as an elongated 
protein stretching away from the cellular membrane with 
HS chains branching out from the unfolded protein back¬ 
bone (1,41,44,49,50). Because MD simulations have been 
remarkably successful in the study of folding processes 
for proteins of similar size (51), we decided to test the ecto¬ 
domain structure using a simulated MD annealing proce¬ 
dure (see Methods). We took advantage of the small 
system size, -9000 atoms, to simulate four HS sequences. 
For each HS sequence, six simulations were performed, us¬ 
ing different initial conformations. The sequence employed 
were: a combination of NA-NA/NS-NS domains (Fig. 2 b, 
Siml-6); a nonsulfonated NA domain (Fig. 2 c, Sim7- 
12); a fully sulfonated NS domain (Fig. 2 d , Siml3-18); 
and a chimeric protonated domain, called NAH (Fig. 2 e, 
Sim 19-24). 

Fig. 3 shows the ectodomain folding for one annealing 
cycle with NA-NA/NS-NS sugar chains. Similar plots 
were observed for all systems and the use of NA, NS, and 
NAH sequences did not lead to any significant difference 
in folding. Fig. 3 a reveals the appearance of secondary 
structure in the ectodomain. There is no protein structure 
during the heating stage, but as the system is cooled down 
and progresses in the 310 K equilibration stage, a-helices 
form at both ends of the ectodomain. 

It is worth noting that the linker positions remain 
unfolded and accessible. Moreover, the structure is still 
highly flexible. Fig. 3 b shows the RMSD calculated from 
a-carbons during the equilibration stage. After 40 ns of 
initial folding, the RMSD oscillates within a wide region 
of 1.5 nm, lacking the plateau characteristic of stable folded 
structures. Our results propose a very mobile ectodomain 
with a combination of folded conformation at both ends 
and unfolded coils at the center, where HS molecules are 
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linked. Fig. 3 c shows a snapshot of the final frame of an¬ 
nealing. Conversely, there was no appearance of compact 
structure for the HS chains, and they remain unfolded for 
all sequence options. 

Fig. 4 shows the average end-to-end distance for HS 
chains during the 310 K equilibration stage. Remarkably, 
such distance is approximately the same for different HS 
sequences—insensitive to the total charge and sulfonation 
pattern. The lack of correlation between sequence and 
end-to-end distance in our annealing simulations is attrib¬ 
uted to three causes: 

1. The HS chains are hydrophilic; therefore, they favor the 
extended conformation in solution. 

2. There is a limitation of the solvent model, the generalized 
Born implicit solvent, which does not compute long-range 
electrostatic forces. Such long-range interaction should 
affect the conformation of a highly charged polymer. 

3. There is a lack-of-crowding effect, in which each anneal¬ 
ing simulation contains only three HS chains. A higher 
HS concentration should force the HS chains to stretch 
due to electrostatic repulsion and steric clashes. 

A quantity of interest is the bending stiffness (El), which 
is widely used in mesoscopic models of glycocalyx (43,52). 
We estimated the El of HS chains using the final stage of the 
annealing cycle, namely the equilibrations at 310 K for 
100 ns. To represent stiffness, we first calculated the persis¬ 
tence length (Pc)- The procedure to compute P c is described 
in the Methods. Then, we employed 


where k B is the Boltzmann constant, and T is 310 K. 

P c was calculated for all four sulfonation patterns. A ta¬ 
ble listing P c and El values is included in the Supporting 



FIGURE 4 Effect of sulfonation pattern on the conformation of HS 
chains. The plot shows end-to-end distance for all HS sequence alternatives. 
Each point represents an average taken over 18 HS chains. Error bars repre¬ 
sent ± SD taken over 18 HS chains. Overlapping regions in the error bars 
highlight the variability in the end-to-end distance values. 


Material. Similar to the end-to-end distances, the P c for 
HS chains was not sequence-dependent. Also, the P c values 
were highly variable, with a mean value of 15.92 nm, in 
good agreement with previously reported P c values for 
all-atom sugar models (37). 

The El mean value for HS chains is 68.28 pN nm 2 , 
approximately an order-of-magnitude lower than 490 pN 
nm 2 , the value commonly used in continuum models. 
Note that to obtain an El of 490 pN nm 2 , the P c value should 
be 114.5 nm, which is unrealistic for HS chains. We 
conclude that the high El value for glycocalyx cannot be 
attributed to the HS chains alone. The other macromolecular 
components of the glycocalyx, such as plasma proteins and 
adsorbed carbohydrate chains, contribute to an enhanced 
stiffness. Although our glycocalyx models include atomic 
information, they still lack the heterogeneity added by the 
adsorbed molecules. 

The second topology block is the transmembrane domain. 
As yet, there is no high-resolution structure available; never¬ 
theless, there is cumulative experimental evidence that all 
syndecans have a-helical transmembrane domains and 
appear as dimers (53,54). In a recent computational study, 
Psachoulia et al. (55) have shown that Syn-2 and Glyco- 
phorin-A share the same transmembrane dimer packing, 
where both a-helices are in contact through a GXXXG 
sequence motif. Because Syn-2 and Syn-4 transmembrane 
domains share 78% amino-acid identity and the GXXXG 
motif (residues 153-157, see Fig. 2 a), we overlapped two 
Syn-4 backbones into a reported Glycophorin-A dimer 
structure. Dimer transmembrane domains were inserted 
into a POPC lipid bilayer. Two systems were built, the first 
with one dimer, the second with three dimers, then inserted 
into a lipid bilayer and equilibrated as described in Table 1 . 

The last topology block, namely the cytoplasmic domain, 
is well characterized experimentally and high-resolution 
dimeric structures are available. During equilibration, the 
structure remained as a dimer, but both ends showed high 
mobility. Certainly, we hypothesize that attachments to the 
transmembrane domain at N-terminus and to actin filaments 
at the C-terminus are necessary to restraint mobility. 

The final stage is to join all three parts into a single 
glycocalyx system, then solvate and ionize it, resulting in 
a system of 5,700,000 atoms. The entire glycocalyx system 
was equilibrated further (Sim29-32). Equilibrium states 
were checked by convergence of temperature, volume, pres¬ 
sure, and total energy. A snapshot of the final frame for the 
three Syn-4 dimer systems with periodic images is presented 
in Fig. 1 c. Movie S1 in the Supporting Material shows the 
NVT equilibration of the three-dimer system. 

Shear-stress effect on glycocalyx conformation 

The endothelial cells are responsive to mechanical stimuli 
from the flowing blood (56), controlling vessel diameter 
(57), vessel permeability (58), and flux of biochemical 
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signals (59,60). The glycocalyx is in direct contact with the 
bloodstream and mediates these stimuli by transmitting the 
flow-induced shear stress and pressure into the interior of 
endothelial cells. In turn, the glycocalyx conformation 
also affects the migration of blood cells and plasma proteins 
through the vasculature. 

We study the effects of shear flow using one of our glyco¬ 
calyx structures, namely the system with three Syn-4 dimers 
and 18 HS chains (Fig. 1 c), which has the highest glycoca¬ 
lyx concentration. To induce shear stress, we apply forces on 
the solvent located in the ectoplasmic side. We performed 
three simulations in the NVT ensemble, using forces of 
0.1 pN (Sim33), 0.01 pN (Sim34), and 0.001 pN (Sim35). 
Such forces result in applied shear stress values of 30.42, 
3.63, and 0.47 MPa, respectively. 

The flow simulations with shear stress of 30.42 and 3.63 
MPa damage the integrity of glycocalyx, as the Syn-4 ecto- 
domains unfold. Fig. S2 shows the final frames of those sim¬ 
ulations. For flow simulation under shear stress of 0.47 MPa, 
Syn-4 ectodomains maintain their structures and only the 
FIS chains are straightened by the flow. Fig. 5 shows the 
initial and final frame as well as the velocity and shear stress 
profiles. The velocity profile is largely parabolic above the 
glycocalyx, but it becomes approximately linear inside the 
sugar layer. Accordingly, the highest shear stress is observed 
at the top border of glycocalyx and it decreases close to the 
membrane. 

We further analyzed this simulation by computing the HS 
charge density and the cytoplasmic domain conformation. 
Fig. 6 a presents the time-evolution of HS charge density. 
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FIGURE 6 Conformational change of glycocalyx for 0.47 MPa shear 
stress flow simulation, (a) The contour plot shows the charge density of 
HS chains in the ectoplasmic compartment. The abscissa represents time; 
the ordinate represents position above the lipid membrane. The values of 
the charge density are color-coded according to the color-bar at the right. 
(b) The plot shows the variation of angle for the three cytoplasmic domains 
of Syn-4 dimers; each line represents one dimer. 
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FIGURE 5 Effect of shear flow on glycocalyx conformation. ( a ) Initial 
conformation of the system: a glycocalyx layer composed of three Syn-4 
dimers {solid) embedded into a lipid bilayer {dark shaded) and 18 HS 
chains {light shaded). The flow is applied on the upper half of the ecto¬ 
plasmic compartment (see text). A graphene barrier was located at 51 nm 
above the lipid membrane (not shown here), {b) Bending of glycocalyx un¬ 
der 0.47 MPa shear stress. {Snapshot) Simulation after 10 ns, periodic im¬ 
ages of the lipid membrane are shown repeated along the x axis, (c) Velocity 
{solid squares) and shear stress profiles {shaded circles). 


Initially, the charge density occupies a region up to 45 nm 
above the lipid membrane. As the flow impacts the HS 
chains, the glycocalyx conformation changes, crowding 
the HS chains into a smaller volume. The charge distribution 
changes accordingly, compressing the negative charges into 
a region of 30 nm within the membrane. Although the flow 
changes the distribution of charges from 45 to 30 nm above 
the membrane, the highest charge concentration remains 
located in a region between 10 and 15 nm above the mem¬ 
brane during the entire simulation (see Fig. 6 a). Further¬ 
more, the simulation reveals that Syn-4 dimers act as 
levers, moving the cytoplasmic domains in the opposite 
direction of the shear flow. To quantify the cytoplasmic 
bending, we measured the angle 8 between each cyto¬ 
plasmic dimer and the z axis. Fig. 6 b shows the cosine of 
8 for the three cytoplasmic domains, where value 1 repre¬ 
sents cytoplasmic alignment with the z axis and value 0 rep¬ 
resents cytoplasmic domain is parallel to the membrane. As 
observed, the three cytoplasmic domains bend toward the 
membrane during the simulation. In endothelial cells, the 
Syn-4 cytoplasmic domains are connected to transduction 
complexes and actin filaments, and it is expected that cyto¬ 
plasmic angling should move the anchoring points in the 
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cytoskeleton, triggering the biochemical machinery inside 
the cell. The Movie S2 shows the flow simulation under 
0.47 MPa shear stress. 

Glycocalyx binding 

Another puzzling aspect of the glycocalyx is its role as a 
flypaper for molecules flowing in the blood; that is, its abil¬ 
ity to recognize and retain extracellular ligands. This is a 
critical function, because ligands deliver biochemical sig¬ 
nals to endothelial cells. An extensive variety of HS-binding 
ligands has been reported, including growth factors, extra¬ 
cellular matrix proteins, cell-cell adhesion receptors, and 
blood-coagulation factors (42). Among the most studied 
ligands are fibroblast growth factors (FGF). The FGF family 
is composed of at least 20 related proteins and it is involved 
in several cellular functions, such as proliferation, migra¬ 
tion, morphogenesis, and angiogenesis (61). In particular, 
the FGF-2 protein is a well-documented case, as several 
high-resolution structures of FGF-2 in contact with HS 
chains are available (62,63). FGF-2 functions by promoting 
dimerization of specific transmembrane proteins, called 
FGF-receptors (FGFR). FGF-2 molecules alone cannot 
gather FGFRs, and HS chains are needed to retain FGF-2 
at the cell surface. As FGF-2s, HS chains, and FGFRs aggre¬ 
gate into macromolecular complexes, FGFRs dimerize and 
transfer the stimulus inside endothelial cells. Nevertheless, 
the molecular mechanism behind the aggregation of FGFs 
and HS chains is still a matter of debate (64). 

We studied the effect of glycocalyx on FGF-2 dynamics. 
In principle, all-atom MD simulations can fully characterize 
the dynamic interplay between ligands and receptors, as it 
has been shown before for systems of up to 100,000 atoms 
(65,66). Unfortunately, the size of our glycocalyx system 
(5,700,000) exceeds that range. Moreover, FGF-2 is more 
massive than typical small ligands, and a thorough equili¬ 
bration between FGF-2 and glycocalyx may require hun¬ 
dreds of nanoseconds or microseconds. To overcome these 
limitations, we employed the Grid-SMD technique, which 
reduces the system size by coarse-graining the atomic de¬ 
tails of glycocalyx into a grid. Then, FGF-2 molecules 
diffuse on the glycocalyx grid; the interacting forces be¬ 
tween FGF2 and glycocalyx are calculated from the grid 
potential. The Grid-SMD approach was initially proposed 
for studies of electric field-driven transport (25), and as 
such it is well suited to study how FGF-2 transport is 
affected by a highly charged polymeric brush. Grid-SMD 
simulations reduce the system size from millions to thou¬ 
sands of atoms, and smooth the energy landscape for fast 
FGF-2 diffusion. The grids that we employed do not take 
into account the steric hindrance of the glycocalyx and 
only consider the mean electrostatic field of the sugar 
chains. Recent developments in the Grid-SMD technique 
allow the addition of extra molecular information. For 
example, steric hindrance can be included by adding a sec¬ 


ond grid potential that represents a three-dimensional poten¬ 
tial of mean force (67,68). 

For binding simulations, we performed an on-and-off MD 
cycle, with the glycocalyx grids acting intermittently on 
FGF-2 molecules under flow. The MD cycle consists of 
three consecutive parts: FGF-2 molecules flowing free; 
then flowing under the effect of glycocalyx grids; and 
then, finally, flowing free again. Such an MD cycle mimics 
the process of protein flowing, binding in the glycocalyx, 
and subsequent release, which can occur by enzymatic 
cleavage of the sugar chains (69). The movie S3 shows 
the on-and-off MD cycle using one-dimer grid. 

Fig. 7 shows the effect of glycocalyx grid on FGF-2’s 
mobility, which is characterized by measuring the FGF-2’s 
average velocity (Fig. 7, a and b ) and average position 
from the membrane surface (Fig. 7, c and d). Fig. 7, a and 
c, refers to one Syn-4 dimer grid (Sim37), whereas Fig. 7, 
b and d, refers to three Syn-4 dimer grids (Sim38). 

For the first part of the cycle, the first 10 ns (Sim36; open 
backgrounds at left), the FGF-2 molecules move in the 
direction of the flow with varying velocities (Fig. 7, a and b), 
and retain their distance from the lipid surface (Fig. 7, c and d). 

The second part of the cycle, the activation of glycocalyx 
grids, is marked by the shaded areas in Fig. 7. FGF-2 mol¬ 
ecules still move in the direction of the flow, but the grid 
forces slow down the proteins and the average velocities 
decrease (Fig. 7, a and b). Remarkably, FGF-2 molecules 
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FIGURE 7 Glycocalyx affects mobility of FGF proteins. The panels 
show flow simulations of FGF-2s turning on and off the electrostatic poten¬ 
tial of glycocalyx. (Open, shaded, and open backgrounds mark states as off, 
on, and off for glycocalyx electrostatic fields, respectively.) (a and b) 
Average velocities for FGFs using one-Syn-4 dimer and three-Syn dimer 
grids, respectively, (c and d) Average positions for FGFs using one-Syn-4 
dimer and three-Syn dimer grids, respectively. Each point represents an 
average taken over 10 FGF-2s. Error bars represent ± SD taken over 10 
FGF-2s. 
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move close to the membrane position (Fig. 7, c and d). The 
FGF-2 molecules concentrate in the region of high HS 
charge density (time zero at Fig. 6 a). These simulations 
suggest that the glycocalyx contributes to gather FGF-2s 
near the lipid membrane, increasing the effective concentra¬ 
tion of FGF-2, which is needed for FGFR activation. 

The third part of the cycle is shown in open backgrounds at 
right. Both simulations were continued, but now the glycoca¬ 
lyx grids were turned off (Sim39 and 40). The FGF-2 veloc¬ 
ities start to increase (Fig. 7, a and by, however, the position 
from the membrane remains unchanged (Fig. 7, c and d). 

It is well established that HS induces not only binding 
but also oligomerization of FGFs; namely, the formation 
of FGF dimers or tetramers along the HS chains (70). 
FGF oligomerization can be an essential step for FGFR ac¬ 
tivity, as it has been hypothesized that FGF oligomers are 
needed for FGFR dimerization (71,72). During our binding 
simulations, we observed the formation of FGF-2 aggre¬ 
gates. Each FGF-2 protein has a charge of +11 e, and con¬ 
tains an amino-acid binding sequence KRTGQYKLGS 
KTGPGQK (each letter refers to one-letter abbreviation 
code for amino acids) to recognize HS chains. As the 
FGF-2 molecules race toward the negatively charged posi¬ 
tions in the glycocalyx grids, they get close to each other 
and assemble into FGF-2 oligomers. To provide a quantita¬ 
tive description of FGF-2 oligomerization, we used an 
iterative nearest-neighbor search procedure to identify ag¬ 
gregates (73). If the centers of masses of FGF-2s are within 
a cutoff distance of 4 nm, the proteins are considered in con¬ 
tact and part of an oligomer. Fig. 8, a and b, shows the FGF- 
2 oligomerization for one- and three-dimer glycocalyx 
grids, respectively. Shaded backgrounds in the x axes 
mark the time region where the glycocalyx grids are active. 
Each circle represents a FGF-2 aggregate; the radius and 
grayscale provide information about the aggregate’s size 
(see legend). Such snowman-like plots provide a schematic 
representation of the aggregation process over time. 

In the first 10 ns, when the glycocalyx grids are not active, 
FGF-2 molecules flow individually across the periodic cell 
and do not oligomerize. When the grids are active (shaded 
backgrounds), monomers crowd near the highly charged 
positions in the grid (see also Fig. 7, c and d), forming olig¬ 
omers of varying sizes. The most remarkable is a short-lived 
pentamer (snapshot in Fig. 8 c), which eventually breaks 
into a dimer and a trimer. For the last 10 ns, the grids are 
turned off again. The newly-formed protein-protein inter¬ 
faces of dimers and trimers remain stable under the flow. 
The existence of similar side-by-side dimers, so-called cis 
dimers, is still under discussion (64), mostly because there 
is no crystallographic evidence yet. It is tempting to state 
that our simulations confirmed the existence of cis dimers. 
However, because the binding simulations are short, the sta¬ 
bility of those oligomers is difficult to assess. We have 
inspected the trajectories of Sim37-40, where FGF2 oligo¬ 
mers appeared, searching for specific FGF2 contacts that 
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FIGURE 8 Glycocalyx induces aggregation of FGF proteins. (a and b) 
Flow simulations of FGFs turning on and off the electrostatic potential of 
glycocalyx. The images show schematic representations of FGF aggrega¬ 
tion under the effect of one dimer (a) and three dimer ( b ) glycocalyx elec¬ 
trostatic fields. ( Circles ) FGF aggregates. For each time frame shown, the 
aggregates are sorted based on their total mass. The number of FGFs in 
each aggregate is coded in the size and grayscale (see inset). (Shaded back¬ 
grounds in the axes) Activation of the glycocalyx electric fields, (c) ( Snap¬ 
shot ) FGF pentamer. (Transparent surface) Isosurface of electrostatic 
potential. 

are known to stabilize in FGF2-FGFR interface (63,74). 
Although we were able to identify the hydrophobic patch 
NYL (N102-Y103-L140 in Plotnikov et al. (74)) as a part 
of the interface in one of the trimers, most of the contact 
surfaces are different. Therefore, it is likely that the olig¬ 
omer contacts are nonspecific and may eventually dissociate 
if simulated longer. 

Altogether, our results provide a generic mechanism of 
FGF-2-HS molecular interplay: The highly charged glyco¬ 
calyx force the intrusion of FGF-2 proteins near the mem¬ 
brane position. Then, the negatively charged positions in 
the glycocalyx act as nucleation centers for FGF-2 oligo¬ 
merization. FGF-2 oligomers are formed within the glyco¬ 
calyx layer. Such aggregates may contribute to trigger 
FGFR activity. Overall, these results show the effect of the 
different grafting densities. The higher the concentration 
of HS chains, the slower the FGF-2s flow and the closer 
they get to the membrane. Also, the high HS concentration 
allow the temporal formation of large aggregates. 

CONCLUSIONS 

We use large-scale MD simulations to present the structure 
and short time dynamics of the glycocalyx layer. We have 
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built the most detailed to-date model of the glycocalyx, 
composed of HS sugar chains connected to transmembrane 
Syn-4 dimers embedded in a POPC lipid bilayer. 

We evaluated four HS sequences containing different 
degrees of sulfonation, including a sequence with regions 
for FGF and antithrombin binding. Our results revealed 
that Syn-4 has a-helical segments in the ectodomain, but 
it is unfolded in the segments where the HS chains are 
linked. Furthermore, HS chains remain extended indepen¬ 
dently of the sulfonation pattern. 

We examined the glycocalyx conformation under shear 
flow. Under moderate shear forces, we observed stretching 
of HS chains, thinning of the glycocalyx thickness, and 
bending of the Syn-4 cytoplasmic domain. This last result 
shows that shear stress can exert influence on intracellular 
structures through mechanotransduction. 

We explored the binding of FGF proteins to the glycoca¬ 
lyx layer. We used Grid-SMD methodology to simplify 
the glycocalyx model, reducing the system’s size from 
5,700,000 to 20,000 atoms. We observed an increase of 
the FGF concentration near negative charges on the HS 
chains, leading to FGF oligomerization. 

Overall, our results provide insights in the structure of 
glycocalyx layer and its dynamics, albeit for nanosecond 
timescales. Our findings motivate the need for further 
improvements; for example, adding different proteoglycans 
and sugar chains, including more atomic complexity into 
the glycocalyx grids, and suitable coarse-graining of 
the glycocalyx network to reach longer timescales. We 
envision that simulations of flow in capillaries that accu¬ 
rately include the effects of the glycocalyx will lead to a 
more comprehensive understanding of blood flow in health 
and disease. 

SUPPORTING MATERIAL 

Directions on obtaining glycocalix structures and Grid-SMD simulations, 
two figures, one table, and three movies are available at http://www. 
biophysj .org/biophysj/supplemental/S0006-3495( 13)01201 -0. 
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